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FOREWORD 



The principal purpose of this paper is to present methods for determinir.g 
minimum-time ship routes on a digital computer. The paper was written in 
two parts. The first treats in detail a numerical method for determining the 
course which requires minimum time to go from one specified place to another 
when the speed of the ship is a known function of position, heading, and time. 

The second part of the paper treats various related topics. The first is 
the corresponding problem in the simpler case wherein the speed is a functic*n 
of position and heading but not of time. The second is a discussion of metho<is 
for obtaining the curves, called equal-time curves or isochrones, which yield 
the maximum distance the ship can attain at any particular time by choosing 
various courses. These are discussed for both the stationary and time- 
varying speed fields; they are rather important because they are easily 
understood and interpreted, and closely parallel numerical methods now in use. 
Third, the problem of effecting rendezvous between two ships is treated. 
Finally a brief discussion is given of various problems and difficulaties which 
may be encountered in computation. 

It has been the intention particularly in the first part to give the numerical 
routines and associated discussions in sufficient detail that a person familiar 
with numerical methods and computers could immediately program and run 
the problem. The method is based on procedures which G. A. Bliss introduced 
in ballistics for calculating differentials. These are applied in a Newton- 
Raphson iteration to determine the course. It is the author's opinion that the 
major outstanding part of the ship- routing problem is the collection of re- 
liable empirical data describing the speed of the ship. 

It is hoped that this paper will serve as an introduction to the adjoint sys- 
tem, calculus of variations, and optimum control theory along lines which aie 
currently being actively pursued in this country and in Russia, including a 
method of solution. The first part of this paper has been accepted for publica- 
tion in NAVIGATION, Journal of the Institute of Navigation. 
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A GENERAL NUMERICAL METHOD FOR DETERMINING OPTIMUM SHIP 

ROUTES 

Frank Do Faulkner* 

A method is given for determining optimum^ ship routes on a 
digital computero This paper is limited to the problem of deter- 
mlni^mlnlmum time„cour,ses^when the speed is a known function of 
the position, heading, and time. 

The method is based on Bliss's methods for calculating dif- 
ferentials, usually applied In a Newton Iteration to determine 
the course. It Is very general, but It requires a knowledge of 
Bliss's adjoint methods. The paper Is selfoontalned; a simple 
case Is given In detail and others outlined. The advent of high- 
speed digital computers In the past few years now makes the solu- 
tion of such problems feasible. 



1. Statement of problem. 

shall take the equations of motion to be expressed in the 

form 

o 

= V cos p 
ly = V sin p 

where x,y are position coordinates, p Is a control variable, 

'' = ''(x,y,p,t) Is a known function of period 211 in p with 
continuous derivatives, t Is time, and the dot C) over a 
variable Indicates Its time derivative. If x,y are rectang- 
ular coordinates, then v Is the speed and p la the angle 
between the velocity vector and t he x axis , the heading angle; 

late Schoor'°'’ "^thematlcs and Mechanics, U. S. Naval Postgrad- 
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If they are, say, latitude and longitude, then there is no such 
simple interpretation. 

We will consider first the problem of going from one specified 
point, (0,0), to another, (X,Y) , with minimum time T. The problem 
is to choose the control variable p as a function of time to ef- 
fect this, and to determine the corresponding curve in x,y,t space. 
In a later section general conditions will be given for minimum 
time courses, as in rendezvous, etc., and corner conditions which 
would be of interest for a sailboat which must tack. 



2. Variational equations, Euler equations. 

In this section some formulas for differentials are derived, 

by procedures which G. A. Bliss introduced in Ballistics. Let us 

consider two neighboring paths, whereon the values of p differ 

r T, , 



by an amount 6p; it is assumed that 



6p| dt is small, 



Then 



the change 5x,6y in x,y satisfy the variational equations 



( 2 ) 



'6x = cos p + Vy6y cos p + (VpCOS p - v sin p)6p 

l6y = sin p + v^by sin p + (v^sin p + v cos p)6p; 



subscripts here Indicate partial derivatives, v = dv/dx, etc. 

Let us multiply these two equations through by two new 
variables (unspecified at present but to be identified a^ 

some stage as Lagrange multipliers), collect terms, and Integrate 
to get 

rT 

(3) [)v(6x - v^&x cos p - Vy&y cos p - Vp6p cos p + v6p sin p) 






= 0 

We may integrate this by parts and rewrite it as 
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(4) [9i6x + y-'5y]Q = 



r T 



[6x(^+Xv^cos p +y^>'v^sln p) + 
+ &y(y^ + XvyCos p +yJVySln p) 



+ 6p(A{v^ ^cos p - V sin p}- +yO'{VpSln p + v cos p}')]dt 

ti+ 



- [v(Xcos p +y^sin p)] 6t^ 



where is a symbol for any point^.or points where p is dis- 

continuous as a function of time, a steering corner . Probably 
there will be no such point for ships and we shall drop the last 
term for the present. To simplify this equation, let us choose 
as solutions to the system 



(5) 



■X +> V cos p + /-^v sin p 

X /X 

* 

+ XVyCOS P + ylv'VySin p 



0 

0 , 



so that the coefficients of 5x,6y in (4) vanish. This defines 
the system (5) which is adjoint to the system (3)» 

In the case of most Interest, the initial values of x,y are 
assumed known; then 6x(0) = 0 = 6y(0) , and (4) becomes 



(6) 0y(T)6x(T) + y(^(T)6y(T) = 

cos p - V sin p) +y»-'(VpSln p + v cos p)]6p dt 

pT- 

= A“v 6p dt, 

Jo P 

-A. ^ 

where A = "Xi and "v = v(l cos p + j sin p). This is the 

fundamental differential formula connecting the change in end 
values of x,y with the variation of the control variable. 

It is a very general condition for ’any optimum that there is 
some solution A to the adjoint system, to be found in the solu- 
tion of the problem, and for that solution, ^ must satisfy the 



condition 
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(7) 



A’v = extreme, 



as a function of the control variable p at every point of the 
course. This condition implies that the coefficient of 6p in 
(6) vanish 



for that solution. Comment. The proof of this condition is es- 
sentially the proof of the fundamental lemma in the calculus of 
variations (see Courant*[l] p 200). Cases wherein it does not 
determine a unique path are easily constructed, but seem to be of 
little practical significance and will not be considered further. 

Equations (5) » (8) are called the Euler equations in calculus 
of variations. The curves whereon equations (1),(5)»(8) are satis 
fled are called extremals of the family (1). Equation (7) or 
(8) may be replaced by the condition that p be chosen so that 



this approach has received considerable attention recently through 
the work of the Russian mathematician Pontrlagin [2] on control. 

The various philosophies of approach represented by (7) » 

(8), (9) are equivalent for this problem, and none eliminates 
the basic problem of solution; namely, conditions are given at 
two values of t, (0,T), with the second one unknown, and a con- 
stant is needed to integrate equations (1),(5)»(8). In the next 
sections a method of solution is given. 

3. Differential formulas. 

Let us consider a fundamental set of solutions for (5) » say 
, *Numbers in brackets refer to references listed at the end 




(9) 




0 



>f the paper 
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taken so that at t = 0 these have the values 



1,0, 0,1, respectively o Then every set of solutions Is a linear 
combination of these two., Now the fundamental equation (6) Is 
linear and homogenepus In the pair so that a multiplicative 

factor will drop out from the formulas for 6x,6y<, Hence we 
may express all solutions of Interest in terms of a single param- 
eter a 



In so far as the use of this formula is concernedo An extremal 

which starts at the origin at t = 0 is characterized by the 

value of a except for the terminal time T, and the problem 

C <3 *^ Ofl^t imuCrvO 

of determining “bourse reduces to that of finding a, To 

It Is shown in the appendix that for extremals a differential 
change in a leads to a differential change 6x(T) ,6y(T) in 
x,y at'- T 



If both a,T are changed the resulting differential changes 
' in the terminal values are 



(10) 





where 



(12) J 




t=T 



( 13 ) 



Ax(T) = x(T)6T - ,^(T)6a 
Ay(T) = y(T)6T + jA(T)6a 
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4. Computational routine « 

Let us guess initially values for a,T, say a^jTQ. 

Then Integrate simultaneously the original equations, a funda- 
mental set solutions to the adjoint system (6), 

with determined then by (10) and p determined by (8); 

calculate also the integral J» 

A curve thus determined will not go to the desired poln^ 

X,Y generally, but to, say, ,Y^ , where i is the itera- 
tion index. In equation (13) , set 
, Ax. = X - X. 

(1^) 

lAy^^ = Y - Y^ 

and determine 6a^,6T^, = T^+bT^^. 

Continue this until some convergence criterion is satis- 
fied, say, 

(X - X^)^ + (Y - Y^)^ < £ , 

where £ is a prescribed number. 

Of course, if a,T are guessed too far from the correct values 
it may not converge; this has not been a problem so far. 

5. General end conditions. 

In the most general case we may have a function g(x,y,t) to 
be minimized, or maximized at the terminal point, subject to N 
constraints of the form 

hj^(x,y,t) = 0, n = 1,‘>'=°N. 

There may be none, one, or two of these (N = 0,1, or 2). In this 
case the cond^jtJ.qns on the en d values , called the transversal con- 
ditions, which must be satisfied for a stationary value of g are 



] (Nx3) 



first (see Bliss [(o] , pZ,0) 

/dh^/dx dh^/dy dh^/dt' 

(15) rank ( dg/dx dg/dy dg/Dt 

\ 'y ^ “)^X 



^ t=T 



and second, the rank of the matrix obtained by omitting either of 
I the last two rows must also be N+lo 

I The first N rows of the above matrix are the coefficients 

of Ax,Ay> 6T in the next row are from Ag, and the last 

row is the corresponding set of coefficients in the differential 
formula (6) rewritten in the form 



(16) 



i'X^x +/^Ay ~(^x-»^y)6T]^^^ 



TT- 

A»v_6p dto 

Jo P 



In the problem studied first, this reduces to the condition 
()ix +yuy)^ 0, and is always satlsfiedo 

One way to determine the path is to guess a,T as before. 
The correctional routines may be obtained from 
Ah^ = Sh^/dxAx + dh^/dyAy + dh^/6t6T, 

and it may be necessary to use 

r A'XiT) = A(T)6T + (dX/da)^6a 
LAyi/(T) = y(^(T)6T + (^/da)^6ao 

The differentials are chosen to drive residual errors to zero. 

Example. Consider the problem of getting to any point 
where x assumes a specified value X, in minimum time, with 
no constraint on y. The matrix of (15) is then 
(10 0 ) 

I 0 0 1 f 

A yO -Ovxy/y)) . 













|r^ 






om 



< 
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Hence yu(T) = 0. This yields one relation 

= (*^^sin a a)6a + yu(T)6T = yu{T) 

for 6a, 6T; the equation for Ax in (13) 
is the other equation neededo 

An alternate procedure is to continue the computation 
until X = X; this determines Then obtained 

by setting ^(T) ^^cos "*”/^2 i®^^ ^i+1 “ seems 

that the only ingenuity required is in guessing starting values 
and in setting up correction routines from (15) » typical comput-» 
ing problems. As this example shows, the correction routines are 
not unique and some may converge better than others. 

6. Comments. 

In a simple example requiring about 100 time steps, each 
Iteration took about one second and a path was obtained in from 
four to twenty seconds on the CDC l604, depending on the function 
V and the initial guesses? this could well be decreased if 
desired. On the other hand, the use of empirical data and the 
attendant calculations will undoubted increase the time. 

There are two other possible methods of solution of this 
problem. If the function v is independent of t, then the 
order of the system may be reduced by eliminating the time t. 

The resulting Euler equation is of order two and may be solved 
by a relaxation procedure (Haltinerj, Hamilton^ 'Arnason [3]),. A 
limited comparison suggests that the differential correction 
procedure converges more rapidly if the number of 

time steps is large. No way is seen to extend the relaxation 
method to differential systems of higher order. The only other 
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method goes by various names: steepest descent or ascent, and 

gradient. It is based on a corrective routine similar to the 
constructive proof of the fundamental lemma of the calculus of 
variations (Courant [l] p 200) » Discussions suggest that all 
have uncertainties associated with convergence. 

Corner condition. If there are any corners (points of the 
curve where v is discontinuous) then the differential formulas 

ti+ 

(6) and ( 15 ) need a term added -[vCXcos p +yusin p)] 6tj^ 

A 

for each point. For the values of A,yu associated with an ex- 
tremal, this vanishes (but does not vanish generally if 
yu =y(^i f etc.). A program is needed to check for this condition, 
if it may occur, since the values of p(t^-), p(t^+) are not 
close to each other. Implicit in the derivation of the above 
correction is that the transients Introduced by "coming about" 
are negligible. 

Many important conditions which must be checked to ensure 
that the resulting curve affords a minimum must be neglected in 
a short paper. These are covered completely , though tersely, 
by Bliss ([6] chapters 7,8). Many of these have no significance 
until a path is obtained. 

As applie d to shin routin g: ^ these have also been discussed 
by de Jong [4P ^contains a discussion of several interest- 
ing cases of particular interest in air navigation. The direc- 
tion perpendicular to /\ defines the curves of constant time 
associated with the extremals from the origin. An alternate 
method of determining minimum courses is to determine these 
isochrones , as done by Hanssen and James [5] 5 , except to use 
extremals and the transversal conditions. This may be done 
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on the computer, as outlined in some preliminary reports 
[lO] and sections 8 and 10 of this paper. 

This method of solution is an application of the methods 
which Bliss [8] introduced for calculating differentials in 
ballistics (see also Bliss [9l> Ch. V for summary, and p. 125 
for other references). 
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Part II OPTIMUM SHIP ROUTING 
Introduction 

Various other aspects of optimum ship routing are taken up here. The 
next sections treat the case where the velocity field does not change signifi- 
cantly with time, and may be considered constant over short periods. The 
differential equations are then of lower order and the problem is simpler. 
Two methods of solution are discussed. Each makes use of extremals, 
which are the minimum- time curves, or brachis trochrones , at least over 
short distances. These may be used together with a relation known as 
'*trans ve rsality” to determine equal-time curves, or isochrones. These de- 
fine the boundary to the points which the ship can reach at any particular 
time: it can generally reach any point up to and on these equal-time curves, 

but no point beyond them at the corresponding time; also points on the iso- 
chrones can be attained only by extremals. The shortest-time route maybe 
determined from these, A second method of solution is given for finding 
the shorte st- time route directly without using the isochrones. 

Next the problem of determining the equal-time curves is taken up for 
the case where the velocity field varies rapidly or significantly with time. 
Then some routines are given for effecting rendezvous between two ships, 

(a) when one ship is following a known course, and (b) when the two ships 
cooperate. Finally a brief discussion is given of problems which may be 
enccuntevel in the computation. 
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/ Constant Velocity Field, Euler Equations, Transversals. 

In this section It will be assumed that the speed v Is 
a function of position and heading, but not time, v = v(x,y,y' 
where x,y are position coordinates, and y' = dy/dx. For 
simplicity, we will generally choose the coordinate set so that 
the Initial point Is the origin and the final point (X,0) 

Is on the x axis. 

The Euler equation, which Is a necessary condition for 
minimum time will be derived now. If x,y are cartesian coord 
Inates, the time required to go a short distance along a curve 
Is approximately 

(17) t = As/v = /l+y'^ x/v , 

and the time to go along any selected curve Is then 



(18) 



T = 



X 



f(x,y ,y' )dx; 



1 2 

this defines f = -\ll+y' /v. If x,y are not cartesian 

coordinates , there are corresponding relations , depending on 
the metric of the coordinate set. 

Let us consider two neighboring paths. On the second, 
y-^y + 6y 



(19) 



y' ■ > y' + 6y' , 



where (6y)' = 6(y' ) , and means "Is replaced by". The 

f I 

Integral |6y'l dx Is "small". The difference In the time 

^0 

required to follow the two courses may be approximated then 

/ X 



( 20 ) 



6T = 



(fy&y + fyi&yMdx + f (x,y ,y' )^^J^6X ; 



subscripts In the Integrand indicate corresponding partial 



-14- 



Pa #32 



derivatives. Let us integrate this by parts to eliminate 6y' 
from the integrand 

( "X X 

(f„ - df , /dx)6y dx + f(x,y,y*)y&X +[f^,6y] . 

0 y y ^ y 0 

If the end points are fixed, then 6X = 6y(0) = 6y(X) =0. If 
T is a minimum, because of the choice of path, then 6T must 
vanish for all allowed functions 6y. Since 6y is arbitrary 
except for being "small" and vanishing at the endpoints of the 
curve, it follows that 



- fy = 0- 



This is the well-known Euler equation and is a necessary 
condition (assuming that f has continuous second derivatives 

We may 

rewrite this as 



and f'yiyi 5^ 0)j see Bliss [6], p. Ilf for a proof 

»/ J 



(23) 



y" = 



fytyi/^y “ ^'^y'y ~ ^y'x-^ * 



which is of the form y" = P(x,y,y'). It may also be rewritten 
as a pair of first-order equations, by introducing a new variable 



(24) 



y' = Z 



z' = F(x,y,z) , 

which is convenient for computations. It should be observed 
that since (23) is of second order, there are two constants of 
integration and an extremal is determined if the values of y,y* 
are given for some value of x; since y(0) = O.ln our problem, 
there is a single constant of integration to be determined. 

The curves defined by (22) or (23) are called extremals . 
Transversals. The equations for curves which correspond 
to equal times along different extremals of a family, such as 
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the family passing through a given point, will be derived now. 
If the endpoint (X,Y) is not fixed, but may vary, then 



the terms in equation (21) associated with it do not necessarily 
vanish. Since we are interested in shortest-time routes, we 
will consider only extremals, so that the integrand and the inte- 
gral vanish. Now the change in the final value of y is 



(25) 



6Y = 6y(X) + y' (X)6X. 



)6y(X) 



6Y 



(see sketch), y' being the value 
on the extremal. Equation (21) 
becomes 

( 26 ) 6T = [f6X -t(6Y - y'6X)fy,]j^ 



J ly'6xy 

6X 

Relations among 
variations of 
end values. 



Now if we have a one-parameter family of extremals such 



as those coming from a given point, then the various points, 
one on each extremal , which correspond to the same value of T 
will be obtained by setting 6T = 0. That is, from (25) and 



(26), 

(27) it - [y’ - f/f , ]jj5X = 0; 

this is the defining relation for the transversal direction 
6Y/6X. A curve S which cuts each of a family ^e]- of extre- 
mals transversally is transversal to the family, or is a 
transversal of the family. Its equation is 



(28) dY/dX = y' - f/fy, ; 

dY/dX is the slope of S and y' is the value on the cor- 

c. 

responding extremal. 

And if we are given a smooth curve S, then there is a 
family of extremals |E| which satisfy (28). The family 
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so defined represents the curves of shortest time from the 
transversal S, at least if the distances are not too great. 



Comment. Sometimes the transversals are perpendicular 
to the extremals. It may be easily verified that this is the 
case for an arbitrary family of extremals if and only if f 



this is the case if v is independent of the heading. 

8. Numerical routine. 

A method of determining a minimum-time ship course is 
now given. It is similar to the method described by Hanssen 
and James [5l (pp259»260) except that it makes use of the 
results Just derived, Involving extremals and transversals. 

Let us take the set of extremals which emanate from the 
starting point. Calculate a family of these with the initial 
heading angle as a parameter, say taking values of tan”^y'(0) 
one degree or ten degrees apart, using the Euler equation. 
Qontinue these for say six or twenty-four hours. The value of 
t must be obtained by integration and the terminal point will 
involve interpolation. At the endpoints determine also, from 
(28) the transveral direction. 

We now have a set of points on the equal-time curve and 
the corresponding direction. Fit a curve, . using this 

data. It is not the most common type of data since not only 
is the point given, but also the slope. Name this curve , 
or case may be. It represents the maximum distance 



has the form 




in the ship-routing problem 
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the ship can be at that time, based on the Initial velocity 
field. 

Now revise or update the velocity field to correspond to 
the values forecast for six or twenty-fours ahead. Then take 
the one-parameter family of extremals which cut , or 

transversally. The slopes of each extremal is given by equation 

(28) , with revised values for v, for the continuing extremals. 
As many extremals as desired can be drawn out from S^. With 
each extremal that is continued from associate the corres- 
ponding value of y'(0). 

Continue these out for another six or twenty-four hours 
and generate another transversal ^48* Update the 

field and continue until a transversal hits the desired terminal 
point. Interpolate to get the initial heading and the route. 

This method has the desirable feature that the results 
may be easily Interpreted, particularly by those already cal- 
culating minimum time routes by present methods. On the other 
hand it requires unnecessary computing to determing one route. 

9. Alternate numerical routine. 

It is felt that the following method will generally 
determine a route more quickly. 

Let us consider equations (24). The equations for the 
variations of an extremal are 
6y' = 6z 

(29) 

6z' = Fy6y + F^6z, 

Let us guess an initial value for y' , Then compute the 
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corresponding extremal , which will probably not go through 

the desired endpoint (X,0). Simultaneously, compute a solution 
to ( 29 ) » with 6y(0) = 0 and any convenient value of 6z(0) , 
since the system is homogeneous and linear in 6y,6z<, The 
value of y(X) will be in error, say it is Treat 6y(X) 

as a difference; set 6y(X) = - and solve for 6z = 6y'(0). 
This gives a corrected value for y'(0). Continue this until 
the terminal value differs from the desired value of y by 
less than some preassigned number, the convergence criterion. 

The values for the velocity field could be updated regularly 
every six hours as before, if desired. A limited number of 
computations suggests that this will give a solution quickly. 

Comment. There is a condition for an extremal to yeld 
a minimum value for the integral, called the non-tangency con- 
dition. If at any point 6y as obtained from (29) above is 
zero, the corresponding point is on the envelope of the family 
of extremals and does not furnish a minimum value to the integral. 
In this case a warning should be given to the operator, or 
a corresponding subroutine generated. This will probably 
never occur for short courses. 

10. Isochrones for time-varying field. 

A method is given here for constructing the Isochrones 
when the velocity field varies with time, and for determining 
the minimum time course by the method of section 8. 

Let us consider various extremals starting from the origin; 
each is determined by a value of a in equation (10) , and 
by equations (1),(5)>(8). At every point of each extremal E 




Interrelations between 
extremal E, velocity v, 
solution A to adjoint, 
and the Isochrone S» 



the relation between the velocity v, the solution A to the ad- 
joint, the extremal E and the Isochrone S is as indicated in 
the figure. Let be the curve defined by the allowed velocity 

V at P, with p as a parameter. 

Then the heading p is chosen so 
that V has a maximum projection 
onto A , see section 2, The 
curve S perpendicular to A at 
P is generated by small changes 
in the parameter a; it is a trans- 
versal of the family. The impor- 
tant property of S is that all 

points on it and to its left in some neighborhood can be reached 
at that time, no points to the left can be, and points of S 
are attained only by extremals. That is, S is part of the 
boundary of a closed region whose points are exactly the points 
where the ship may be at that time t. 

We then get points on a curve S and the tangent direction. 
These define the isochrone S(t), A relatively small number of 
points is required to define S since the tangent direction is 
also given; this tends to be offset by the fact that v is given 
from empirical data. If t is small enough, the curve S is 
similar to a circle or an arc thereof. If the isochrone S is 
given, it must be emphasized that the normal to it determines 
A , and that the extremals are not generally perpendicular to S, 

The first Isochrone S(t) which touches the specified terminal 
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point obviously yields the minimum-time route. The value of 
T and the parameter a must be determined by Interpolation, 

If the total time exceeds the time that forecasting can be 
done reliably, then the extremals may be replaced by terminal 
segments of great circle routes, the extremals for a uniform 
velocity field, as done by Hanssen and James [ 5 ] (p 26I) , or 
statistical data may be used as the basis for determining v. 
Whatever data is "most reliable" should be used; the results 
can be no more reliable than the data, 

11. Rendezvous between ships. 

Two ships will be said to rendezvous if at some time 
their positions coincide. It is assumed that the time required 
for the terminal maneuvering is negligible compared to the 
total time spent. 

In the first example, let the position of the second 
ship be denoted by x*,y*, these being known functions of 
time. As before (see section 4), we guess values for the 
parameters of the extremal, a,T. In the correction routine 
for these we must allow for the distance the second ship 
will travel if we change T by an amount 6T. In place of 
equations (13),(l4) , use 

X* - X = (5; - x*)6T - 

( 29 ) . 1 

y* “ y = (y - y*)6T +')\J5a, 

for corrections 6a, 6T to the estimates of a,T, all quan- 
tities in (29) being evaluated at t = T after the first 
estimate of the trajectory has been calculated, J was defined 
by equation (12). 
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The corrected values for a,T are used and another 
course computed. This Is continued until the distance be- 
tween the two is less than some preassigned number, at the 
end of a computation. 

Cooperating ships. The problem is more interesting 
if the second ship cooperates. Let quantities associated 
with the second ship be denoted by an asterisk (*). For 
rendezvous, there must be a time T such that x(T) = x*(T), 
y(T) = y*(T). There is a further condition that 

(30) A(T) II A*(T) ; 

the vectors defined by the adjoint systems must be parallel 
at time T. That is, if H = then 

(31) H = 0. 

We will need also the differential change in in 

our computational routine. If we change a,T by small amounts, 
we get a differential change in the terminal values 

f 

aA = [-Aj^(T)sln a +7^2 (T)cos a]6a + A(T)6T 

(32) 

= [•^j^(T)sln a aj6a + y>'(T)6T 

and two more , exactly like these, with starred terms, asso- 
ciated with the second ship. 



Computational routine. Guess a,a*,T and calculate a 
first approximation to the courses for the two ships. Let 
the values at the end of the I'th iteration be x^^ ,y^^ ,x* ,y* . 
Generally x^^ x* , or y^^ y* , or A is not parallel to A" 
at t = T, Use the differential formulas as differences to 
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drive the errors or residuals to zero: 



(33) 



= (Xj^ - x^)6T - J^ba + J^*6a* 
y* - = (yj^ " y*)6T + j'Xba - J*>,*6a^ 



-(y* 

+ [(-A]^sln a a^* 

-[(-P\*sin a* +;>^2COS a*J^ 



- pM -y^X*)6T 

-(70'lSin a yj^cos a)A*]&a 

-(y^*sln a* yJ^Qos a*)'\ ]6a*. 



If this converges, it yields a stationary value of T. The 
computation would be continued until some convergence criterion 
is satisfied, say until (x^-x*)^ + (y^^-y*)^ + < £ , where 

c is a preassigned number. 

An alternate procedure is the following. Suppose the co- 
orlnates are chosen so that the first ship is initially at (0,0) 
and the second at (X,0), with X > 0. Guess a, a* as before 
and compute the two extremals , stopping when at some time t = 
T^, X = X*. Equations (33) again constitute a set of three 
equations for three unknowns, though only two of these, 6a, 6a*, 
are needed to start the next iteration. 

An alternate method would be to determine the Isochrones 
S(t) ,S*(t) for each ship, by the method described in section 
10, continuing until they touch. Since (or if) they are smooth 
curves , S and S* are tangent to one another at the first 
time of contact. This is equivalent to condition (3D» that 
AHA*, since A and S are alvjays perpendicular. 
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12. Comments 

As remarked initially the first purpose of this paper has been 
to give a method for determining optimum ship routes.. It may have 
given however an impression that the problem is simpler than it 
actually is. In this section a brief discussion is given of some 
situations which may arise, which are not treated earlier. 

It should be remarked that the routine given converges to 
a stationary value of time, which may not be a minimum. This is 
particularly the case for long routes, where a course not near 
that chosen may yield a smaller value of time. There is some con- 
dition which is an extension of the envelope condition in elemen- 
tary calculus of variations, and the well-known conditions of 
Weierstrass and Legendre (see Bliss [6], chapter 8). This does 
not lessen the value of this paper, since these conditions can 
be checked only after a course has been determined. 

The envelope condition in the simplest case states that if 
an extremal furnishes an extreme value to an integral, then it 
cannot contain a focal point or a point of the envelope of the 
family of extremals from the initial point. For the simplest 
problems this can be checked by checking to see if 6y is zero 
at any point of the course, but no simple check is known generally 
to the author. 

If the course is determined by the method of isochrones, then 
it is suspected that the Isochrone will then develope an interior 
corner. Unfortunately the method of fitting the Isochrones will 
probably have an implicit assumption of smoothness which may well 
obscure the existence of the corner. The difficulty is aggravated 
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by the fact that the data is given empirically and we must try 
to smooth out or minimize the random errors by a curve-fitting 
or smoothing routine. 

For the initial part of the course we will surely be able 
to predict accurately. For the later part, we may be required to 
use statistical data, or in the absence of that, a terminal sec- 
tion consisting of part of a great circle. If the times involved 
are large, so that better data becomes available, the course 
should be continually redetermined, using the current position 
as t1ie starting point. This is particularly important on long 
routes, and when storms Introduce large changes in the speed 
which cannot be estimated accurately beforehand. 

It seems that time may well be too simple a cost function. 
Cargo damage, danger, passenger discomfort, etc., increase with 
wave height. The cost function must be made up by someone 
familiar with the details of shipping. In the best routing 
procedures, several activities must be coordinated. (1) First, 
the ^meteorologist or oceanographer must collect the data and 
predict sea state, and weather, too, if it is significant. 

(2) A cost engineer familiar with the ship and its cargo must 
make up a realistic cost function in terms of time, danger, 
etc. (3) The mathematical framework is given here. (4) Finally, 
there is no little problem of programming properly. Some study 
needd to be made of methods of smoothing the data, so that the 
derivatives, which are notoriously poor when the data is poor, 
are satisfactorily smooth. The success of the method depends 
on the Integrated efforts of all of these. 
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Finally, let us consider a problem where the above routine 



breaks down. Let us consider a sailboat In a narrow channel 
where hills on either side slow down the wind. Let us consider 



and w are constants. It Is easily verified that the maximum 
speed upwind (In the direction of the positive x axis Is v^/4. 
However, to effect this, the boat must come about an Infinite 
number of times; this Is known In control theory as the "chattering" 
solution. The above-given routine is too simple for this problem: 
the time spent coming about and the ensuing transients are not 
negligible. It seems to the author that some feel for the 
mechanics of such problems is necessary in the programmer, that 
one must suspect beforehand that such events are likely. These 
seem to be of little importance with powered ships and were not 
investigated in any detail. 
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speed In the form v = v^(l - y /w )(1 - 2 cos p) , where v^ 
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Appendix. In this section the differential formulas 
of section 3 are derived. Let Ai,A 2 denote as before the so- 
lutions to the adjoint system which have initial values i,j, 
resp. Let a = a^cos a + AgSin a, as before and a = |a|. 

Then = Aj^cos^a + Ag^sin^a + 2 A^'Agsin a cos a. Let 
X = Vcos Pj y = V sin p, and denote the "velocity" by V; 
lower case \^v • was used earlier. Lower case symbols will be 

—4 “♦ -4 

used to denote unit vectors in this section, V = Vv , a = 

AW, Then 

7 = v^v + Vn, 

p p » 

-4 — ^ —4 

where n = kxv. If u is the unit vector parallel to V^, 
then 

u = k X a/a 



and - 

V = /(V 2+V2) u 

-4 -4 

on every extremal, since V^^a s 0. 

Now let us consider at the same time and place two ex- 
tremals differihg by small amounts 6p,6ao Then 



-4 

(A,l) ^pp*'^ 

Now also 



So 

(A.2) 





IVpIpU IVpIkxu 

iVpIpk!^- !Vyiw> 

= -/(V2+Vp2) A. 



^4 

^ II «AW 



Now consider also 



^The terms to the right of the parallel bars || in the 
right-hand margin will be used to suggest the operations by 
which the next equation is obtained. 
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V »«A = 0 = V »(aiCOs a + Aosin a) 
P P ^ 



Let also 


P* = 


(Vp-Al)®+ (Vp-Xj)2 j 


then 


sin a 


= -Vp*Ai/p , oos a = Vp-Aj/p 


and 






(A.3) 


V .A 
P a 


= Vp»(-Aisin a + Agcos a) 



= Po 



Hence from (Ad) ,(A,2) ,(A„3) 

6p = p6a/(A/[V2+Vp2]), 

and the differential formula (6) derived earlier becomes 
when A = Ai 

»T 

(Xi6x + Ui6y)m = - (p^sin a/A)dt da 

1 Jq 

for extremals. But 



p2 = [y(V2+Vp2)^-Aj^ + L/(V2+Vp2)2^»A2]* 

= ^ {[(Ajcos a + Agsin a) "Aixjtjz 

a2 '■ 

+ ['(a^cos a + Agsin a)«A 2 xIc] 2 | 



Hence 



and 



V2+V 2 



Xi Xg 
dl U2 



I 

(Xj6x + ^^6y)^ = -J 



(X 26 X + M 26 y)^ = J 



-J 


pT y2^y 2 


Xi X 2 




'0 A^ 


Ml M2 


= J 


pT y2+y^2 


Xi X 2 




'0 A^ 


Ml M2 



dt sin a 6a , 



l|-d2 



dt cos a 6a , 
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Hence 



6x(T) = - J w(T) 6a 
6y(T) = J X(T) 6a, 

the values of X,u being those associated with the maximizing 
condition and 



J = 



^2 

Ui Ug 1 



V2+V 2 


Xi, Xg 


‘'O 





dt 



This derivation seems long and complex but no simplification 

has been found, and the corresponding three-dimensional formula 

has so far defied derivation. ^ 

Frank D, Faulkner 6 June 1962 
Monterey, Calif. 
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